Strand diffusion-limited closure of denaturation bubbles in DNA 



o 



o 



I 

o 
o 



(N 

q 

en 

o 

> 

X 

J3 



Anil Kumar Dasanna, Nicolas Destainville, John Palmeri, and Manoel Manghj^] 

(1) Universite de Toulouse; UPS; Laboratoire de Physique Theorique (IRSAMC); F-31062 Toulouse, France, EU 

(2) CNRS; LPT (IRSAMC); F-31062 Toulouse, France, EU 

The closure dynamics of a pre-equilibrated DNA denaturation bubble is studied using both Brow- 
nian dynamics simulations and an analytical approach. The numerical model consists of two semi- 
flexible interacting single strands (ssDNA) and a bending modulus which depends on the base-pair 
state, with double-strand DNA (dsDNA) segments being 50 times stiffer than ssDNA ones. For 
DNA lengths from Ai" = 20 to 100 base-pairs (bp) and initial bubble sizes of A — 6 bp, long closure 
times of 0.1 to 4 /xs are found, following a scaling law in A^'*. The bubble starts to close by a fast 
zipping which stops when the bubble reaches a highly bent metastable state of length around 10 bp. 
The limiting final step to complete closure is controlled by the dsDNA "arms" rotational diffusion, 
with closure occurring once they are nearly aligned. The central role of chain bending, which cannot 
be accounted for in one-dimensional models, is thus illuminated. 

PACS numbers: 87.15.H-,87.15.A-,82.39.Pj 



I. INTRODUCTION 

Understanding the dynamics of biological processes 
such as transcription, duplication or DNA translocation 
by viruses is a challenge for biophysicists. By itself, the 
dynamics of double-stranded DNA (dsDNA) is a complex 
issue which has been tackled at two different scales : (i) at 
the macromolecule scale, dsDNA is a model semi-flexible 
polymer whose internal structure is ignored and whose 
dynamics is controlled by thermal fluctuations and bend- 
ing modes [I1I2]; (ii) at the base-pair scale, the dynam- 
ics focuses either on the base-pair closure during DNA 
renaturation at room temperature, the so-called DNA 
zipping |3], or on the breathing dynamics, i.e. the fast 
opening and subsequent closure of small bubbles, a lo- 
cal opening of successive base-pairs, with very low prob- 
ability [IHZ]- Such studies do not consider the closure 
dynamics of a thermalized or pre-equilibrated denatura- 
tion bubble which couples both scales, i.e. local base- 
pair closure and chain diffusion modes. The closure of 
such a denaturation bubble occurs for instance at the fi- 
nal stage of DNA transcription, when RNA-polymerase 
leaves the locally open DNA |5], or is observed in vitro 
in DNA solutions as a rare event with the largest time 
scale. Altan-Bonnet et al. measured the closure dynam- 
ics of denaturation bubbles of 18 base-pairs (bp) by fluo- 
rescence correlation spectroscopy and found surprisingly 
long closure time in the 20 — 100 /iS range [9] . Such long 
times were also measured in bulk experiments on hairpin 
formation in ssDNA and RNA oligomers, which are much 
slower than theoretically estimated times of end-to-end 
contacts for simple semi-flexible polymers [10]. Bubble 
lifetimes of about 1 /iS have also been observed for DNA 
oligomers of 14 bp in NMR measurements of the imino 
proton exchange [TT| . 

Several models have been used for studying bubble 
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breathing and in attempts to explain these large ex- 
perimental bubble lifetimes, (i) The Poland-Scheraga 
model |13j is a one-dimensional (ID) Ising model mod- 
ified to account for the entropic penalty of creating 
a closed fiexible loop. This term leads to a non- 
monotonic free energy landscape in which the typical 
breathing time comes from a Kramer's process [At. i5ul2]. 
(ii) The Peyrard-Bishop model is a non-linear phonon 
ID model where bubbles emerge as soliton-like solu- 
tions of undamped Newton's equations in a Morse poten- 
tial [71 dH [T3] . Although these models capture the short 
time scale breathing dynamics, they are not adapted 
to the issues of renaturation and "equilibrated" large 
bubble closure since they miss DNA diffusion in solu- 
tion. But strand dynamics is expected to be dominant 
at least for long DNAs and large bubbles, since for a 
flexible chain of length N the Rouse diffusion time scales 
like A^ ^ . (iii) Other models fit well the experimental auto- 
correlation function [21 HSl [H] , but with the relaxation 
time as a fitting parameter, which does not shed light on 
the origin of such large times. 

In an attempt to explain the mechanism behind these 
large bubble lifetimes, we focus on the out-of-equilibrium 
closure of a thermalized denaturation bubble using both 
Brownian Dynamics (BD) simulations and analytical ap- 
proaches. We implement two numerical models where 
the different bending rigidities of dsDNA segments with 
a persistence length of roughly £ds = 150 bp and ssDNA 
ones with £ss = 3 bp are explicitly included, and whose 
coupling with base pairing has been central to under- 
standing equilibrium properties [TH [19] . We show that 
the denaturation bubble closure occurs in two steps (Fig- 
ure u_ I . The first step consists in a fast zipping of the 
initia bubble until a metastable bubble state of length 
^ 10 bp is reached. The driving force for this fast ki- 
netics is the energetic gain in base-pair closure at room 
temperature. At some point, the metastable bubble is so 
bent that zipping becomes forbidden by the large bending 
energy cost to close the bubble. The second step of the 
closure is then controlled by the relaxation of the bent 
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FIG. 1: Snapshots of a typical Brownian dynamics simulation 
showing (a) the initial equilibrated bubble, (b) the metastable 
state, and (c) the bubble just before closure (A'^ — 40 bp). 
(d) Sketch of the metastable state. 



state through the diffusion of the two stiff dsDNA arms. 
Once the two rigid strands are aUgned, bubble closure 
occurs almost instantaneously. The first mean passage 
closure time is found to scale with the DNA length N (in 
bp units) as Tdosure ~ Ar2-4±o.i foj. 20 < N < 100 and 
mesoscopic parameters typical of real DNA. 



II. MODELS 

We simulated the closure of a large bubble of initial 
length Lq = A'^ — 6 in the middle of an homopolymer 
DNA. We used both BD simulations of two interact- 
ing semi-flexible strands, and the Kinetic Monte Carlo 
(KMC) algorithm which simulates the mean semi-flexible 
chain with an internal Ising spin dynamics corresponding 
to the bp state (broken or unbroken) [T8[ [19] . 



A. Brownian dynamics simulations 

The DNA is modeled by two interacting bead-spring 
chains each made of N beads located at r, . The Hamil- 
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where tj — rj+i — Vj and ij = tj/|tj|. The first term of 



the rhs. of Eq. (IT]) is the stretching energy with stretch- 
ing modulus f3Ks — 100 {f3~^ — k^T where T is room 
temperature) and a ~ 0.34 nm is the equilibrium dis- 
tance between two beads in each strand. The second 
term is the usual bending energy with a bending mod- 
ulus K(,j that depends on the local chain configuration 
(£p = /Sntj). The interaction energy between the two 



strands (the Hydrogen bonding between two complemen- 
tary bases) is modeled via a Morse potential p4j of width 
A and depth A: 
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where pj = jr^ — r^- | is the distance between comple- 
mentary bases at position j along the chain and po — 
1 nm is the equilibrium distance. The stacking interac- 
tion is modeled by a bending modulus Kf, which depends 
on p, interpolating from Kds/2 — 75 fee T for dsDNA state 
to Kss — Sk^T for single stranded one, according to [6] 
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where f{pj) = [1 - 
transition and pb 



erf(^^)]/2. A' is the width of the 



1.5po- The variable bending mod- 
ulus depends on three consecutive base-pair distances, 
which provides cooperativity. We chose A = 0.2 nm, 
A' = 0.15 nm, po = I nm, and f3A = 8.^ The thresh- 
old value for p, discriminating between open and closed 
states, is fixed at 1.13 nm (a slightly different value does 
not change the results). 

The evolution of r , (t) is governed by the overdamped 
Langevin equation 



C^ = -^rM{rj}) + Ut) 



(4) 



where C = Sirria is the friction coefficient for each bead of 
diameter a {rj — 10""^ Pa.s is the water viscosity), ^j(i) is 
the random force (with zero mean), which mimics the ac- 
tion of the thermal heat bath and obeys the fluctuation- 
dissipation relation (^j(t) • ij{t')) — &k^T(^5ij5{t — t'). 
The adimensional time step, 5t — Stk-^T / {a'^ i^) , was 
fixed, for sufficient accuracy, at 5 x 10^** {St — 0.045 ps). 
The initial bubble is created by turning off Hint and then 
equilibrated for 2 ps. Output values are then calculated 
every 10^ steps once Hint is turned on, and total simula- 
tion times range between lO"^ to 10^ steps (0.4 to 4/is). 
Samples are made of 200 trajectories and error bars are 
standard deviations. 



B. Kinetic Monte Carlo simulations 

Out-of-equilibrium dynamics of the mean DNA chain 
(the center of mass of the two strands) has also been 
explored numerically by KMC simulations. We imple- 
mented the coupled model defined in Refs. [IHlIin] where 
the mean chain is composed of N identical beads rep- 
resenting the base-pairs. Simulation details are given 
in Ref. [5D] (each bead now represents one base-pair 
and has the mobility of a pair of beads in BD simula- 
tions). At each Monte Carlo Sweep of physical duration 
St = 0.019 ps, a bead is chosen at random and a random 
move is attempted for this bead. In addition, at each 



■li I I iiiri"^'"^**^"iiirr 




FIG. 2: Top: melting map (position of open base-pairs in 
DNA as a function of time). Bottom: time evolution of the 
bubble length, L{t)/Lo, in red, the bubble end-to-end dis- 
tance, R{t)/{aL(t)), in purple, and the tangent-tangent prod- 
uct C{t) in blue (A'^ — 60, time is in units of 10'* BD steps 
= 0.4 ns). 



Monte Carlo step, we also attempt to flip the sign of one 
Ising spin variable Ci, according to a standard Metropolis 
procedure. However, it might be that, in a real DNA, the 
frequency of change of internal degrees of freedom is dif- 
ferent from this arbitrarily chosen one. To rule out this 
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possibility, we simulated various systems where 10 
100 spin-flips are attempted per St. The average closure 
times then changed by at most ±20% as compared to 
1 spin-flip, thus proving that this is not a critical issue. 
The good matching between the KMC and BD results 
below supports this observation. 



III. CLOSURE DYNAMICS 

The evolution of the bubble size (normalized by the ini- 
tial bubble size Lq), L{t)/Lo, is plotted as a function of 
time in Figure [2] for a typical simulation {N = 60). Two 
regimes appear clearly: at short times, L{t) decreases 
rapidly which corresponds to a zipping dynamics, un- 
til it reaches a metastable state characterized by a small 
bubble of constant size L ~ 10, the center of which dif- 
fuses slowly (cf. Fig. [2] top). The simulation is stopped 
when the bubble closes, which defines the Mean First 
Passage Time (MFPT) for closure. The dimensionless 
bubble end-to-end distance R{t) / aL{t) and the tangent- 
tangent product C{t) = tiniftend, where ti„it (resp. tend) 
is the mean value of the two strand tangent vectors at the 
beginning (resp. end) of the bubble (see Fig.Il]), are also 
plotted. First, R{t) < aL{t) during the entire simulation 
run, except at closure when R{t) ~ aL(t) which favors 
the formation of dsDNA which is stiff at this length scale. 
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FIG. 3: Log-log plot of the total closure time, Tdosurc, vs total 
DNA length N. Dots corresponds to BD simulations results 
and triangles to KMC ones. The solid line is a fit between 
A = 50 and 100 yielding Tdosure ~ A^"*^. 



Likewise, C{t) starts here from a positive value and then 
decreases to a negative value in the metastable state close 
to —1, which corresponds to anti-correlated tangent vec- 
tors. It then undergoes large fluctuations, and the bub- 
ble closure corresponds to a value of C(Tciosuro) ~ 1, i-e. 
when the two stiff arms are aligned. This type of be- 
haviour is observed whatever the initial condition C(0). 
The zipping is faster than the diffusion time of the small 
arms so that the distance between both DNA extremities 
remains almost constant while the bubble "pushes" in the 
direction parallel to the arms such that C(rzip) ~ —1. 
The DNA adopts an "hairpin" configuration as observed 
in Fig. [Tb. The correlation between C{t) and L{t) is a 
clear indication that the spatial configuration of DNA 
plays a central role in bubble closure. 

Closure times simulated using both BD and KMC are 
shown in Figure [3J For one, the two types of simula- 
tions yield almost the same closure times and variation 
with N, showing that this quantity does not depend on 
the specific model. Since in KMC simulations the two 
strands are not simulated explicitly, this indicates that 
closure dynamics is dominated by the whole chain dy- 
namics, which is the slowest process. For another, the 
fit of numerical data yields Tdosuro ~ ^" with an ap- 
parent exponent a = 2.4 ±0.1 for the largest values of 
N and almost a = 3 for small N values. In any case 
an exponent a > 2 is found which is larger than a = 1 
found in breathing dynamics using the Poland Scheraga 
model [H |S], a = 0.52 in thermal renaturation [3], or 
a = 1.37 in anomalous zipping dynamics [5T]. We argue 
below that this apparent exponent is the signature of a 
complex dynamics governed by the rotational friction of 
the (almost) rigid arms in the metastable state. We now 
discuss in details the two regimes. 
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FIG. 4: Log-log plot of the zipping time (circles) and the 
metastable (MS) residence time (triangles) vs the zipped ss- 
DNA length P = 3{Lo — L)/5 and the arm's length M respec- 
tively. Fits lead to r^ip ~ P^* and Tmot ~ Af^ to M^'". 



A. Fast zipping process 

In Figure [4] it is shown that the fast zipping process 
is characterized by an anomalous exponent for L{t) and 
the zipping time 
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required to zip P = |(^o ~ L) bp. We also simulated 
DNA renaturation with an open end (data not shown) 
and found a zipping regime alone, with the same expo- 
nent. A similar "anomalous" exponent (a ~ 1.37) has 
been found in a numerical study of the zipping dynamics 
of two polymers connected at one end [21j . This anoma- 
lous dynamics is interpreted, by analogy with polymer 
translocation [211 [22] , as an out-of-equilibrium propaga- 
tion of the tensile force located at bubble ends, / ~ A/a, 
along the bubble backbone (the analog of the force ap- 
plied on the monomer located at the tiny pore in polymer 
translocation) [23]. The dynamics of L{t) is governed 
by the equilibrium between friction, driving and bending 
forces along the chain 
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where h{t) is the bubble segment size involved in the fric- 
tion, and i^bend is the bending energy stored in the bub- 
ble. Since the initial bubble is pre-equilibrated, -Ebcnd 
can be neglected and the initial zipping is governed by 
the two first terms of Eq. JgI. For a fully flexible bubble, 
it has been proposed [23l - l25 l that b{t) ~ aL{tY (where 
V is the Flory exponent) which yields r^jp ~ P^+'^/J, 
an exponent value slightly larger than 1.4. Note that 
(i) as the tensile force straightens the bubble, the bending 
force increases thus defining a moderate forcing regime. 
In this regime (a/i?o < Paf < 1), Sakaue found, for 
the translocation case P^, an exponent value of 1.43 



which agrees well with Eq. ((sl). (ii) Moreover, at the 
end of zipping, the ssDNA is not fully flexible since 
L « 2£ss, which slightly increases the friction and thus 
the anomalous exponent. Finally, the bubble geometry 
imposes that the zipping process stop when the bubble is 
highly bent. From Eq. J6]) the metastable state is reached 
when £'bond(^) — fo-L. For a (circular) bent bubble, 
Ehcnd{L) = 2Tr^2Kss/L which leads to Z ~ ^/Airhi^ijA, 
i.e. a few base-pairs. In other words, when L{t) — L, 
the bending energy cost for closing one more bp becomes 
larger than the base-pairing gain. 



B. Diffusion limited closure of the metastable state 

The metastable residence time Tmot is plotted in Fig- 
ure |4] as a function of the arm size, M. Clearly a scaling 
law Tmet — Af^ appears for the longest arms. For shorter 
(and thus stiffer arms), the exponent is larger, around 
2.4. This law might be surprising since /^k^s/M ~ 3 — 10 
and one would expect the arms to be stiff. A scaling esti- 
mate for two stiff arms connected by a flexible joint yields 
a time inversely proportional to the rotational diffusion 
coefficient of a rodlike polymer of length M (neglecting 
prefactors and logarithmic terms coming from hydrody- 
namics) [21] 
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To check this law numerically and since we are limited in 
M due to computational cost, we performed BD simula- 
tions with a larger dsDNA persistent length, /Sk^s — 400. 
Figure [5] shows that both the metastable residence time 
and the closure time now scale like r^ct — M^ and 
Tciosuro ~ N^ respectively. This thus demonstrates that 
the closure process is limited by the rotational diffusion 
of the arms. The reason why the exponent deviates from 
3 to 2 when we consider a real DNA with (Sk^s = 150 
is twofold. First, the arms being semi- flexible, bending 
modes enter into play which decrease slightly the friction 
and accelerate the dynamics. Second, the rotational dif- 
fusion law Eq. (l7| is only valid at long time scales and the 
closure may occur before the diffusive regime is reached. 
To clarify this behaviour, we model the system as fol- 
lows. The dynamics of closure is characterized by two 
coordinates, namely the end-to-end distance of the two 
single strands forming the bubble of length L, R{t), and 
the angle between the tangent vectors located at bubble 
ends, 9{t) = arccosC(i) (see Fig.ll]). The metastable res- 
idence time, Tmet(r) [where r = {R,6)], is defined as the 
MFPT needed to go from this metastable state to the 
state where R ~ 2aL and ~ 0c ^ 1- It is solution of a 
backward Smoluchowski diffusion equation in a potential 

U{r)m 

- (3S/U ■ V[DT„,et(r)] + A[i?T,„et(r)] = -1 (8) 

and scales like r^^t ~ D^^ and T,„ot ^ R^/D'^ for ro- 
tational and translational diffusion, respectively. Solv- 
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FIG. 5: Log-log plot of (i) the metastable residence time and 
the closure time vs the arm's length M (circles) and chain 
length A'^ (triangles) respectively for /?Kds = 400. Fits lead to 
Tmet — M^ and Tdosurc ~ N^ ■ (h) The metastable residence 
time vs M for /3Kds = 30 (squares) showing a saturation for 
M > 2/3 Kds. 



ing Eq. ([8]) is out of reach, and we first concentrate 
only on rotational diffusion. The potential U{9) = 
Kss(l — cos9)/L, which favors the 6 = state, is the 
effective bending free energy of the bubble. We find 
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where Ei[z] = 
flexible bubble. 



- J S^dt. In the limit of an infinitely 
L/dss -^ CO, Eq. ([9]) simplifies into 
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Simulation data show that, in the metastable state, the 
probability distribution of 6q has a maximum around 
2.5 rad. By choosing a small value such as Oc = tt/IO as 
observed and L ~ 24s, one finds Tmct — t,^j/2. Eq. ^ 
reproduces qualitatively the decrease of Tmot when Kgg 
increases, as observed in Fig. |6] for /3kss = 1,3 and 7. 
However, the inset shows that the argument leading to 
Eq. ([9]) does not reproduces quantitatively this variation. 

The diffusion of the bubble should also be taken into 
account. We discuss two limiting cases: (i) a stiff bent 
bubble which corresponds to a maximal correlation be- 
tween bubble and arm rotational diffusion. For a circu- 
lar bubble, one has geometrically R = 2aLsm{9/2)/9 
and the metastable residence time is still given by 
Eq. ^. (ii) The opposite, fully flexible bubble (4s < L) 
where the two diffusion processes are uncorrelated. The 
metastable residence time is roughly the longest time of 
(i) T^^t defined in Eq. (fTol) and (ii) the MFPT necessary 



FIG. 6: Closure time probability for /3kss = 1,3 and 7 (300, 
600 and 400 samples respectively) and A^ = 100. The average 
value, Tmct, is plotted in the inset and clearly decreases with 
increasing Kbb- The solid line corresponds to Eq. Kt. 



to go from the initial size R{Tzip) to an almost extended 
configuration R ~ aL, which for an entropic spring with 
I3U{R) = R^I{Na^), is the Kramers-like result [S^: 



T]f3a^L{4:L + M) 



dx / dy 
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where both bubble and arm translational frictions come 
into play. As soon as L ~ 10 < M this MFPT is 
smaller than r^^ . In the real case of a semi-flexible bub- 
ble (Z ~ ^ss), the two diffusive processes are entangled, 
which yields larger closure times in N"^'"^, as observed in 
the simulations. 

Finally, although simulations were done for stiff arms 
with M < /?Kds, we have checked by decreasing /3Kds 
from 75 to 30, that for long semi- flexible arms, T,„ot(Af) 
saturates for M > 2/3 k^s, as shown in Fig. [5J Indeed, to 
close the bubble it is sufficient that the two segments of 
length ~ ^ds adjacent to the bubble align since the rest 
of the dsDNA chain is decorrelated. This yields an upper 
limit for the diffusion time of t^^^ ~ 7//?(2a^ds)^ — 20 /us 
for 4s = 150 bp. 



IV. DISCUSSION 

Using Brownian dynamics simulations of a coarse- 
grained model of DNA, we show that the closure of a 
pre-equilibrated denaturation bubble of size larger than 
10 bp is dominated by the rotational diffusion of the stiff 
dsDNA arms. The first stage consists in a fast zipping, 
the driving force of which is the base pairing energy on 
the order of SksT (TS], leading to an average growth 
rate of 10^ bp/s (even though the time dependence of 
the bubble size is not linear). This fast zipping regime 
stops when the bending energy stored in the small bubble 
is measured to be roughly 20 k^T. The closure is then 



difFusion-controUed with a small driving force due to the 
bending torque applied by the bubble. We numerically 
find closure times Tdosuro ^ A^^* for DNA lengths smaller 
than the dsDNA persistence length of 150 bp. For larger 
DNA lengths, we show that the closure time saturates at 
20 fj,s. 

Our key conclusion is that DNA closure dynamics is 
governed by chain bending, which is a 3 dimensional ef- 
fect that cannot be accounted for by ID models (such as 
Poland-Scheraga [4j [Sj [13] or Peyrard-Bishop [H] mod- 
els). Aside from small time scales, much smaller than the 
Rouse time of the DNA chain (where DNA breathing dy- 
namics is governed by the opening/closure of base pairs) , 
DNA bubble dynamics involves several intertwined in- 
gredients such as base-pair closure, DNA chain diffusion, 
and elastic forces. 

Taking into account hydrodynamic interactions would 
slightly decrease the apparent exponent of the scaling 
laws. For instance, the rotational diffusion time of a 
rigid rod, given in Eq. ([7|, would be divided by a term 
in In M [26] . Hence as soon as small and thus stiff DNAs 
are considered, the correction to the free draining case re- 
mains small. We did not consider the helical structure of 
DNA in this work. In recent numerical work |28J, Baiesi 
et al. show that the unwinding dynamics of a helical 
polymer made of two interwound strands (without any 
attractive interaction) leads to a unwinding time which 



scales as iV^-^^. One might expect that the closure of a 
denaturation bubble would include both winding of the 
bubble plus rotational diffusion of the two arms in or- 
der to relax the superhelical stress [IS]- The precise role 
of the axial rotational dynamics of the two interwound 
strands is currently under study. A lower bound for the 
closure time can be estimated by considering the winding 
of an initial bubble of length Lq. The stochastic evolu- 
tion of a winding angle (/)(i) = 2TrL{t)/p (where the DNA 
pitch is p = 10.5 bp) is given by Crot(Od0/di = n{t) 
where n{t) is a random torque along the DNA axis with 
(n(i)n(i')) = 2kBT(:,ot{t)Sit - t'). Since Crot - i(t)"" 
where zv — 2, |, | (Rouse, Zimm models in theta and 
good solvent respectively [26]), one finds after integra- 
tion 
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2+zv 
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where 3.5 < 2 + zv < A. Hence, DNA winding might be 
the missing mechanism explaining the quantitative gap 
between our estimate of the closure time of 0.04 /is for 



a 30 bp long DNA (see Fig. 
measured experimentally [5] 



and values of 1 — 100 /is 
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